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Abstract 

A new 8-node decohesion element with mixed- 
mode capability is proposed and demonstrated. The 
element is used at the interface between solid finite 
elements to model the initiation and propagation of 
delamination. A single displacement-based damage 
parameter is used in a strain softening law to track the 
damage state of the interface. The method can be used 
in conjunction with conventional material degradation 
procedures to account for inplane and intra-laminar 
damage modes. The accuracy of the predictions is 
evaluated in single mode delamination tests, in the 
mixed-mode bending test, and in a structural 
configuration consisting of the debonding of a stiffener 
flange from its skin. 

Introduction 

Delamination is one of the predominant forms of 
failure in laminated composites due to the lack of 
reinforcement in the thickness direction. Delamination 
as a result of impact or a manufacturing defect can 
cause a significant reduction in the compressive load- 
carrying capacity of a structure. The stress gradients 
that occur near geometric discontinuities such as ply 
drop-offs, stiffener terminations and flanges (Figure 1), 
bonded and bolted joints, and access holes promote 
delamination initiation, trigger intraply damage 
mechanisms, and cause a significant loss of structural 
integrity. 

The fracture process of high performance 
composite laminates is quite complex, involving not 
only interlaminar damage (delamination), but also 
intralaminar damage mechanisms (e.g. matrix cracking, 
fiber fracture). For effective predictive capabilities, 
progressive failure analysis tools for all modes of failure 
are needed. 



Figure 1. Experiment illustrating stiffener- flange 
debonding. 

The objective of this paper is to present a method 
to simulate progressive debonding or delamination 
based on decohesion elements. A criterion for mixed- 
mode delamination is proposed. The effectiveness of 
the method is assessed with an emphasis on its ease of 
use and the accuracy of the predictions. 

Numerical Simulation of Delamination 

The study of delamination mechanics may be 
divided into the study of delamination initiation and the 
analysis of delamination propagation. Delamination 
initiation analyses are usually based on stresses and use 
criteria such as the quadratic interaction of the 
interlaminar stresses in conjunction with a characteristic 
distance 1, 2 . The characteristic distance is an averaging 
length that is a function of geometry and material 
properties, so its determination always requires 
extensive testing. 

Most analyses of delamination growth apply a 
fracture mechanics approach and evaluate strain energy 
release rates G for self-similar delamination growth. 
The G values are usually evaluated using the virtual 
crack closure technique (VCCT) proposed by Rybicki 
and Kanninen 3 . The VCCT technique is based on 
Irwin’s assumption that when a crack extends by a small 
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amount, the energy released in the process is equal to 
the work required to close the crack to its original 
length. The Mode I, Mode II, and Mode III energy 
release rates can then be computed from the nodal 
forces and displacements obtained from the solution of 
a finite element model. The approach is computationally 
effective since the energy release rates can be obtained 
from only one analysis. 

In the present paper, an approach is proposed that 
is well suited to nonlinear progressive failure analyses 
where both ply damage and delaminations are present. 
The approach consists of placing interfacial decohesion 
elements between composite layers. A decohesion 
failure criterion that combines aspects of strength-based 
analysis and Fracture Mechanics is used to simulate 
debonding by softening the element. The proposed 
constitutive equations for the interface are 
phenomenological mechanical relations between the 
tractions and interfacial separations. With increasing 
interfacial separation, the tractions across the interface 
reach a maximum, decrease, and vanish when complete 
decohesion occurs. The work of normal and tangential 
separation can be related to the critical values of energy 
release rates. 

Decohesion can be implemented as a material 
response such as the Dudgale-Barenblatt (D-B) type 
cohesive zone 4 . The D-B cohesive zone model was first 
applied to the analysis of concrete cracking by 
Hillerborg et al 5 . The concept has also been used by 
Needleman 6 to simulate fast crack growth in brittle 
solids. Needleman considered that cohesive zone 
models are particularly attractive when interfacial 
strengths are relatively weak when compared with the 
adjoining material, as is the case in composite 
laminates. 

In order to predict the initiation and growth of 
delaminations, an 8-node decohesion element shown in 
Fig. 2 was developed and implemented in the ABAQUS 
finite element code. The development of this element is 
based on prior work 7, 8 . The decohesion element is used 
to model the interface between sublaminates or between 
two bonded components. The element consists of a 
zero-thickness volumetric element in which the 
interpolating shape functions for the top and bottom 
faces are compatible with the kinematics of the elements 
that are being connected to it. The material response 
built into the element represents damage using a 
cohesive zone ahead of crack tip to predict delamination 
growth. The concept of interface elements has been 
used in different types of problems: compression after 
impact 7 ' 9 , damage growth from discontinuous plies 10 , 
and diametrical compression of composite cylinders 11 . 


Constitutive Equations 

The need for an appropriate constitutive equation 
in the formulation of the interface element is 
fundamental for an accurate simulation of the 
interlaminar cracking process. A constitutive equation is 
used to relate the stress <7 to the relative displacement 8 
at the interface. Some strain softening models that have 
been proposed are shown in Fig. 3 and include: linear 
elastic-perfectly plastic; linear elastic-linear softening; 
linear elastic-progressive softening; linear elastic- 
regressive softening; and Needleman 6 . One 
characteristic of all softening models is that the 
cohesive zone can still transfer load after the onset of 
damage (<f in Figure 3). For pure Mode I, II or III 
loading, after the interfacial normal or shear stresses 
attain their respective interlaminar tensile or shear 
strengths, the stiffnesses are gradually reduced to zero. 
The area under the stress-relative displacement curves is 
the respective (Mode I, II or III) fracture energy. Using 
the definition of the J integral proposed by Rice 12 , it can 
be shown that for small cohesive zones, 

a( 5)d5 = G c (1) 

where G ( is the critical energy release rate for a par- 
ticular mode, and S F is the corresponding relative dis- 
placement at failure (8 PP , S pro , Si in or 5 re in Figure 3). 



Figure 3. Strain softening constitutive models. 


Bilinear Softening Model 

The linear elastic-linear softening (bilinear) model 
is the simplest to implement, and is most commonly 
used 7 9, n ’ I3 ' 15 . The double cantilever beam test shown 
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Figure 2. Eight-node decohesion element; t~0. 


in Figure 4 illustrates the material response. Point 1 in 
Figure 4 is subjected to a low tensile load that is within 
the linear elastic range. A high initial stiffness K p 
(penalty) holds the top and bottom faces of the interface 
element together. Point 2 represents the onset of 
damage. In single-mode delamination, the stress at point 

2 is equal to the corresponding interlaminar strength of 
the material, cr r . As the relative displacement increases, 
the interface accumulates damage and the stress is lower 
than the strength (point 3). The energy released at point 

3 is the area of the triangle 0-2-3. The energy released is 
the most important component of the failure criterion 
described in the next section. If the load were to 
reverse, point 3 would unload to the origin, as shown in 
the figure. 

The critical value of the energy release rate is 
attained at point 4. For any relative displacement larger 
than point 4, the interface does not carry any tensile or 
shear loads (point 5). In other words, at point 4 all the 
available interfacial fracture energy has been 
completely consumed. Note that when modeling 
delamination with a softening response, the 
delamination tip is not defined explicitly. While the 
onset of damage occurs at point 2 in Figure 4, the 
delamination tip could be defined as the point where the 
tractions at the interface are zero, which is point 4. 

The softening response illustrated in Figure 4 is 
representative of the tension or the shear response but 
not compression. It is assumed that compression loads 
do not cause delamination or softening, and the effect of 
compression on damage of the interface was neglected 
in the present work. 



Figure 4. Bilinear constitutive model. 


The concept of decohesion zones to simulate 

delamination growth in composites is usually 

implemented by means of interface elements connecting 
the individual plies of a composite laminate. 

Decohesion elements can model the discontinuity 

introduced by the growth of delaminations. They can be 
divided into two main groups: continuous interface 
elements and point interface elements. Several types of 
continuous interface elements have been proposed, 
ranging from plane interface elements with zero 
thickness connecting solid elements 7 ' 9 ; plane interface 
elements with finite thickness connecting shell 
elements 11 ; line interface elements 10, [5, 16 ; and spring 
interface elements that connect pairs of nodes’ 3, 14 . 

The bilinear interfacial constitutive response shown 
in Figure 4 can be implemented as follows: 

i) 8 <8° => the constitutive equation is given by: 

<r = K p 8 (2) 

ii) 8° <8 < 8 F => the constitutive equation is given 
by: 

a = (l-D)K p S (3) 

where D represents the damage accumulated at the 
interface, which is zero initially, and reaches 1. 
when the material is fully damaged. 

iii) 8 > S F => all the penalty stiffnesses is set equal to 
zero. If crack closure is detected, interpenetration is 
prevented by reapplying only the normal stiffness. 
Frictional effects are neglected. 

The properties required to define the bilinear 

interfacial softening behavior are the initial stiffness 

(penalty) K P , the fracture energies G /c , G//o and Gwc 
and the corresponding nominal interlaminar tensile and 
shear strengths, T and S. The accuracy of the analysis 
depends on the penalty stiffness K P that is chosen. High 
values of K P avoid interpenetration of the crack faces 
but can lead to numerical problems. Several values have 
been proposed for the penalty stiffness, K P : 10 7 N/mm 3 
[Ref. 9], 5.7x 10 7 N/mm 3 [Ref. 16], 10 8 N/mm 3 [Ref. 
13]. Other authors have determined the value for the 
penalty stiffness as a function of the interface 
properties. Daudeville et al. 17 have modeled the 
interface as a resin rich zone of small thickness, t h and 
proposed a penalty stiffnesses defined as: 


The problem of contact of the crack faces after 
failure is addressed by re-applying the normal penalty 
stiffness. The process of reapplying the normal stiffness 
when interpenetration is detected is typical of solution 
procedures of contact problems using penalty methods 
in a constrained variational formulation. 



t, 


w 

p 



(4) 


where G 2 3 , Gj it and E 3 are the elastic moduli of the 
resin rich zone. After a substantial number of numerical 
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experiments, Moura determined that a penalty 
stiffnesses of only 10 6 N/mm 3 for all modes produces 
essentially the same results while avoiding potential 
convergence problems during the nonlinear procedure. 

In order to fully define the interfacial behavior, the 
unloading response must be specified. Petrossian et al. 16 
have proposed an unloading curve with a slope 
corresponding to Hookes law. Such a procedure, 
typically used in the formulation of plasticity problems, 
would lead to the use of the same penalty stiffness when 
reloading and to permanent relative displacements 
along the interface when the load reverts to zero. 
Crisfield et a!. 15, 16 and Daudeville 17 , on the other hand, 
have proposed that with reversing strains the material 
unloads directly toward the origin, as shown in Figure 
4. The assumption is that during reloading the 
interfacial stiffness is lower than the original 
(undamaged) stiffness. Such a procedure simulates the 
effects of the previous damage mechanisms that 
occurred along the interface and was therefore adopted 
in the present work. 


The delamination mechanisms in Mode II and 
Mode III are assumed to be the same. Therefore, Mode 
III can be combined with Mode II by using a total tan- 
gential displacement 5p defined as the norm of the two 
orthogonal tangential relative displacements 5 X and Sy as 

s n =fif7sj ( 6 ) 

The total mixed-mode relative displacement S m is 
defined as 

S m =j8f+Sf, (7) 

where S z is the relative opening (Mode I) displacement. 
Using the same penalty stiffness in Modes I and II, the 
interlaminar stresses are 

a z =K P S l 

r xz = K P S x (8) 

T yz = Kp Sy 


Mixed Mode Delamination Criterion 

In structural applications of composites, 

delamination growth is likely to occur under mixed- 
mode loading. Therefore, a general formulation for 
interface elements must deal with mixed-mode 
delamination growth problems. 

Under pure Mode I, II or III loading, the onset of 
damage at the interface can be determined simply by 
comparing the stress components with their respective 
allowables. Note that the onset of damage does not 
imply the initiation of delamination, since the tractions 
closing the crack at onset are at their maximum value. 
Under mixed-mode loading, however, damage onset 
may occur before any of the stress components involved 
reach their respective allowables. 

A mixed-mode criterion is proposed here that is 
based on a few simplifying assumptions. Firstly, it is 
assumed that delamination initiation can be predicted 
using the quadratic failure criterion 



where o z is the transverse normal tensile stress and t xz 
and Ty Z are the transverse shear stresses. T and S are the 
nominal normal tensile and shear strengths, respec- 
tively. 


The single-mode failure initiation displacements are 
then 


Sf =T / K P 
8l=S!K P 


( 9 ) 


where T and S are the nominal tensile and shear 
strengths of the interface. If the relative opening 
displacement S z is not zero, the mode mixity can be 
expressed by 


( 10 ) 

S l 

The mixed-mode damage initiation displacement is 
obtained by substituting Eqs. 6-10 into 5, which gives: 


c0 cO 
d i °a. 


\ + P 2 


(<?,°) 2 + ( M u ) 


0x2 


(ID 


A quadratic interaction between the energy release 
rates was selected. The interaction law defines the total 
interfacial fracture corresponding to the point where the 
interface is unable to transfer tensile or shear loads. The 
quadratic interaction criterion can be expressed as 



( 12 ) 
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where the equality is only satisfied at failure. For the 
bilinear stress-displacement softening law assumed 
here, the critical energy release rates in Mode I and 
Mode II are simply the areas under the triangle 0-2-4 in 
Figure 4. 


G/r = 




(13) 


The Mode I and Mode II energies released are 
computed from the area of the shaded triangle 0-2-3 in 
Figure 4, which is 


G, =G* 


G// - G uc 


6?S?K p {S t -S?) 



<*?- 

)K P \ 

si) 

(<V 

-<;) 

1 


-si 

) 


(14) 


where Sf and S[ f are the ultimate opening and tangen- 
tial displacements, respectively. The ultimate opening 
displacements are calculated using Eqs. 9 and 13. Using 
the definition of the mixed-mode relative displacement 
5 m in Eq. 7 and the mode mixity ratio P given by Eq. 10, 
one can solve for the ultimate relative displacement 

by substituting Eqs. 14 into the equality of Eq. 12. After 
some manipulation, the result is a quadratic equation 
that has only one positive root, which is: 


S F J 1 + Z x 

S r *S r M F > a -(*2) 2 + 2 fit# - JW + \ 


( 15 ) 


triangle 0-C-D is the bilinear strain response in Mode I. 
It can be observed that the tensile strength T is lower 
than the shear strength S , and the ultimate displacement 
in shear is lower than in tension. In this three- 
dimensional map, any point on the 0-1-2 plane 
represents a mixed-mode relative displacement. 

Stress 



Figure 5. Combined plot of single-mode bilinear 
material responses. 


The map of all softening responses under mixed 
mode is illustrated in Figure 6. The curve FI represents 
the stresses resulting from the displacements at the 
onset of damage given by Eq. 11, while the curve 
labeled G represents the ultimate relative displacements 
calculated with Eq. 15. The triangle 0-A-B is the 
bilinear softening law for a mixed-mode relative 
displacement of S m . The triangle 0-A-B is identical to 
the triangle 0-2-3 in Figure 4. For reference, the triangle 
0-C-D in Figure 6 is the Mode I bilinear softening 
response. It can also be observed that the effect of 
compression on the material response is neglected. 


where S? F =8 f -S? and sf, F =S[, -Sf n which are 
calculated using Eqs. 9 and 13 and the following 
required material parameters: the penalty stiffness K p , 
the interlaminar strengths T and 5, and the material 
toughnesses G {c and Gj jc . 

In summary, the mixed mode softening law pre- 
sented above is a single-variable response similar to the 
bilinear single-mode law illustrated in Figure 4. Only 
one state variable, the relative displacement variable 5 m , 
is used to track the damage at the interface. By 
recording the highest value attained by S m , the 
unloading response is such as shown in Figure 4. The 
displacements for initiation 5 Q m (P) and ultimate failure 
8 F m {p ) are functions of the mode mixity P and are 
computed with Eqs. 1 1 and 15, respectively. 

A mixed-mode softening law can be illustrated in a 
single 3D map by representing Mode I on the 2-3 plane, 
and Mode II in the 1-3 plane, as shown in Figure 5. The 


Stress 



Figure 6. Map of strain softening response for mixed 
mode delamination. 
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Element Formulation 


The element stiffness matrix is based on the 
standard isoparametric linear Lagrangian interpolation 
functions for three-dimensional (8-node) elements. The 
relative displacements between the top and the bottom 
faces of the element in a local coordinate frame x-y-z are 



= BU 


(16) 


where B is the matrix relating the element’s degrees of 
freedom U to the relative displacements between the top 
and bottom interfaces. 

The three-dimensional form of Eq. 3 is 


integral over the area of the element, which gives the 

^ 7 o 

following element stiffness ’ : 

K,, fm =j> r ((l-D)C)B<M (21) 

The integration is performed numerically using a 
Newton-Cotes integration, which has been shown to 
perform better than Gaussian integration in problems 
involving strain softening 13, l4 . The integration points of 
a zero-thickness decohesion element coincide with the 
four corners of the element. Since the material softening 
response is evaluated at each integration point, the 
element can soften one corner at a time, giving it the 
potential to model non self-similar delamination growth. 


Nonlinear Solution 



a. 


a 

o = (i -D)C8 or • 

T xz 

II 

1 

O' 

n 

S x 


.V 


Sy 


(17) 


where I is the identity matrix, C is the undamaged 
constitutive matrix 


C = 


K P 

0 

0 


0 0 
K p 0 
0 K p 


(18) 


and D is a diagonal matrix representing the damage 
accumulated at the interface: 


D = 


d 0 
0 d 
0 0 


0 

0 

d 


(19) 


The term d on the diagonal is the damage parameter, 
which is a nonlinear function of , the highest 
mixed-mode relative displacement experienced by the 
material 


rf = 


8 F m {* 

'max c0 | 

m ~°m ) 

o max 

- *m) 


( 20 ) 


Using the maximum value of the relative 
displacement rather than the current value prevents 
healing of the interface. 8™ ax is the only state variable 
that needs to be stored in the database to track the 
accumulation of damage. 

The minimization of the potential energy subjected 
to the kinematic constraints of Eq. 17 leads to the usual 


The nonlinear solution of the problems presented 
here was performed using standard ABAQUS 
procedures. However, the softening nature of the 
interface element constitutive equation causes 
convergence difficulties in the solution of the analysis. 
The solution of most problems requires tens of load 
increments averaging 200-300 iterations per increment. 
Such a large number of iterations could render the 
solution of most common problems computationally 
impractical. Schellekens 14 recently suggested that in 
problems where failure is highly localized the 
displacement norm in Riks method should be 
determined considering only the dominant degrees of 
freedom. May 18 describes a new automated solution 
procedure for structures with strain-softening materials 
that is based on a constraint equation that uses only the 
displacement parameters associated with the localized 
failure zone in such structures. Unfortunately, a local 
arc-length procedure was not available for the analyses 
presented here. 


Results and Discussion 

Four test problems were selected to validate the 
decohesion elements. The first problem consists of the 
double cantilever beam (DCB) test used to determine 
Mode I toughness. The second problem modeled 
consists of the end-notched flexure test (ENF) used for 
Mode II toughness. The third test is the mixed bending 
mode test (MMB). All three of these problems have 
analytical solutions that were developed by Mi and 
Crisfield 19 . These closed form solutions provide an 
approximate framework against which to assess the FE 
models. The fourth test case examined is a structural 
problem of skin and stiffener flange debonding. 
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DCB Test for Mode I 


The ASTM standard specimen used to determine 
the interlaminar fracture toughness in Mode I ( G lc ) is 
the double-cantilever beam (DCB) specimen. This 
specimen is made of a unidirectional fiber-reinforced 
laminate containing a thin insert at the mid-plane near 
the loaded end. A 15-cm.-long specimen, 2 cm.-wide, 
and composed of two 1.98-mm-thick plies of 
unidirectional material was tested by Morais 20 . The 
initial crack length is 5.5 cm. The properties of the 
graphite material are shown in Table 1, and the 
properties of the interface are shown in Table 2. 


Table 1. Properties of the Graphite/Epoxy material 20 . 


En 

E22=E33 Vi2=Vl3 V23 

Gl2=Gl3 

G 23 

150. MPa 

11. MPa 0.25 0.45 

6.0 MPa 

3.7 MPa 

Table 2. 

Properties of the DCB specimen interface 20 . 

Gic 

Gnc T 

S 

K p 

0.268 N/mm 1 .45 N/mm 30. MPa 

40. MPa 

1 0 6 N/mm 3 


Using Eqs. 9 and the properties of the interface 
shown in Table 2, the relative Mode I and Mode II dis- 
placements for damage onset are t) 0 =3 0.x 10 6 mm. and 
<5° =40.xl0' 6 mm., respectively. The corresponding 
ultimate relative displacements calculated from Eqs. 1 1 
are Sf =17.9x10 3 mm. and Sf t =72.5xl0' 3 mm. 

The ABAQUS finite element model, which is 
shown deformed in Figure 7, consists of two layers of 
C3D8I incompatible-mode 8-noded elements. C3D8I 
elements are superior in bending to other low-order 
continuum elements. The anticlastic effects were 
neglected and only one element was used across the 
width. One hundred and twenty elements were used 
along the span of the model shown in Figure 7. 

A plot of reaction force as a function of the applied 
end displacement d is shown in Figure 8. The beam so- 
lution was developed by Mi and Crisfield 19 for isotropic 
adherend materials and using plane stress assumptions. 
Note that the beam solution is somewhat stiffer than the 
test and FEM results which is probably due to the 
assumption of isotropy in the analytical solution. After 
the initiation of delamination, fiber bridging in the test 
specimen causes a small drift in the response compared 
to the FEM and analytical solutions. 



Numerical studies with different element sizes 
indicate that the accuracy of the prediction can be 
significantly lower if the size of the elements used in the 
softening zone is greater than some maximum value. 
The maximum predicted load sustained by the DCB 
specimen calculated using several mesh densities is 
shown in Figure 9. The results indicate that poor results 
are obtained for this problem when the element size is 
greater than 1.25 mm. This mesh size is consistent with 
the results of Gonsalves 9 , who used l-mm.-long 18- 
node quadratic elements for the analysis of a DCB 
specimen. 



Figure 8. Load-deflection response of DCB test. 
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120 r 


300 r 



20 - 

0 1 1 t 1 J 

0 0.5 1 1.5 2 2.5 

Element Length, mm. 

Figure 9. Debond load as a function of element size. 
ENF Test for Mode II 

Even though the end-notched flexure test specimen 
shown in Figure 10 exhibits unstable crack propagation 
for short crack lengths, its simplicity makes it a 
common test to measure Mode II fracture toughness. 
The length of the specimen modeled here is 10 cm., its 
width is 1 cm., and the initial debond length is 3 cm. 
Aluminum adherends were used rather than composite 
to achieve a closer approximation to the analytical 
solutions calculated by Mi 19 . The thickness of the 
adherends is 1.5 mm. The properties of the interface are 
the same as for the DCB model. 

. Applied load 



Figure 10. ENF test specimen. 

The load-deflection responses for the finite element 
model and the analytical prediction are shown in Figure 
11. It can be observed that both solutions are in 
excellent agreement. 

Mixed Mode Bending Test 

The most widely used specimen for mixed-mode 
fracture is the mixed-mode bending (MMB) specimen 
shown in Figure 12, which was proposed by Reeder and 
Crews 21 and later re-designed to minimize geometric 
nonlinearities 22 . The main advantages of the MMB test 
method are the possibility of using virtually the same 
specimen configuration as for Mode I tests and varying 
the mixed mode ratio from pure Mode I to pure Mode 
II. 



Figure II. Analytical and FEM load-deflection curves 
for ENF test. 



Figure 12. MMB test specimen. 

The specimen analyzed here has a length of 10 cm., 
a width of 1 cm, and an initial debond length of 3 cm. 
The thickness of the aluminum adherends is 1.5 mm. 
The properties of the interface are the same for the DCB 
and ENF models. The length of the MMB lever e was 
chosen as 43.72mm, which corresponds to a ratio of 1 
for G/G fI and to a ratio of 2.14 between the load at the 
mid-span of the beam and the opening load. The MMB 
load fixture is simulated by applying an opening load of 
100 N. at the edge of the debond, and an opposite load 
of 214 N. at the mid-span of the beam. 

The model is composed of two layers of 100 
C3D8I solid elements. A deformed plot of the finite 
element model is shown in Figure 13. 



Figure 13. Deformed plot of MMB Model. 
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The load deflection curve calculated from FEM and 
two beam solutions from Ref. 19 are shown in Figure 
14. One beam solution assumes a quadratic interaction 
between the energy release rates, and the other assumes 
a linear interaction. No test results were available for 
the configuration analyzed. Reasonably good 
correlation is obtained among all three methods. 



Displacement, mm. 

Figure 14. Load-deflection plot for MMB test. 

Skin-Flame Debonding 

As was mentioned earlier in this paper, the 
objective of the present work is to develop a 
methodology that can predict non-self-similar 
delamination growth without user intervention such as 
for remeshing or for releasing constraint equations. The 
effectiveness of the methodology for modeling common 
structural problems without unusually high 
computational or modeling costs is assessed. 

To evaluate the use of decohesion elements in 
structural problems, a problem involving skin-stiffener 
flange debonding was selected. Extensive testing of this 
problem has been performed by Cvitkovich and 
O’Brien 23, 24 , and analyses have been conducted by 
Krueger and O’Brien 25 . The configuration of the 
specimens studied in Refs. 23 and 25 is shown in Figure 
15. The specimens are 203 mm. -long, 25.4 mm. -wide. 
Both skin and flange were made from IM6/3501-6 
graphite/epoxy prepreg tape with a nominal ply 
thickness of 0.188 mm. The skin lay-up consisting of 14 
plies was [0/45/90/-45/45/-45/0] s and the flange lay-up 
consisting of 10 plies was [45/90/-45/0/90] s . 



Reference 25 describes the complexities of this 
problem, such as the presence of several interacting 
failure modes, a non self-similar crack growth, crack 
fronts that jump between plies and do not remain at a 
single interface. References 23 and 25 also show that 
delamination seldom initiates or grows between the 
flange and the skin, but rather on a surface one or two 
plies away from the skin-flange interface. A photo 
micrograph taken from Ref. 23 and shown in Figure 16 
shows these features. 

The properties of the unidirectional graphite/epoxy 
and the properties of the interface reported in Ref. 25 
are shown in Tables 4 and 5, respectively. 



Figure 16. Detail of debond area in a failed specimen 24 . 
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Table 4. Material Properties 25 . 


IM6/3501-6 Unidirectional Graphite/Epoxy. 


En 

E22=E33 

Vi2=V]3 

V23 

Gi2=Gi3 

G 23 

144.7 MPa 

9.65 MPa 

0.3 

0.45 

5.2 MPa 

3.4 MPa 


Table 5. Properties of the interface 25 . 


G,c 

Guc 

T 

S 

0.075 N/mm 

0.60 N/mm 

61 MPa 

68 MPa 


The relative Mode I and Mode II displacements for 
damage onset calculated using Eqs. 9 and the interface 

properties shown in Table 5 are <5,° =61.xl0' 6 mm. and 
=68.xl0‘ 6 mm., respectively. The corresponding 
ultimate relative displacements calculated from Eqs. 1 1 
are Sf =2.46x1 O’ 3 mm. and = 17.5x10 3 mm. 

A complete analysis of the delamination growth in 
the tension specimen requires a high degree of 
complexity. For instance, Krueger 25 developed highly 
detailed two-dimensional models using up to four 
elements per ply thickness. Krueger modeled discrete 
matrix cracks and delaminations including the specific 
paths followed by the delamination in particular section 
planes. 

The approach taken here, on the other hand was to 
determine if it is possible to predict the debond load of 
the specimen using decohesion elements in a much 
coarser three-dimensional model. To keep the modeling 
difficulties low and the approach applicable to larger 
problems, the model that was developed uses only two 
brick elements through the thickness of the skin, and 
another two through the flange. It is implied from such a 
simplified geometry that the interaction between the 
skin and the flange is the physical mechanism that leads 
to the debonding of the flange. The complete model 
consists of 1,002 three-dimensional 8-noded C3DI 
elements and 15,212 degrees of freedom. To prevent 
delamination from occurring at both ends of the flange 
simultaneously, model symmetry was reduced by 
modeling the tapered end of the flange with a refined 
mesh on one side and a coarser mesh on the other. This 
model does not contain any pre-existing delaminations. 

Unsymmetric layered properties were assigned to 
the three-dimensional elements by using the composite 
solid section option in ABAQUS. The option calculates 
the elemenfs stiffness by performing a Simpson’s 
integration through the thickness using material points 
in each of the plies within the element. Gaussian 
integration is used in the lamination planes of the 


element. Approximations to the lamination planes were 
necessary in the tapered sections of the flange because 
planes of integration cannot be defined parallel to the 
plies in the tapered elements. 

Deformed plots of the finite element model 
immediately before and after flange separation are 
shown in Figure 17. It can be observed that only the 
refined end of the flange separates. A detailed analysis 
of the results indicated that the coarse end of the flange 
also softens, but that the separation of the flange at the 
refined end relieves the stresses at the coarse end. It is 
interesting to note that static tests exhibit debond of 
only one end. Fatigue tests, on the other hand, induce 
debonding at both ends. 

Note that the debond growth is not symmetric: the 
debond initiates on the left corner of the flange shown 
on the top of Figure 17 due to the unsymmetry 
introduced by the terminated plies at the flange tapered 
ends. Cvitkovich also commented on the difference that 
he observed in the modes of failure of the two corners 23 . 



Figure 17. Deformed plots of skin/flange debond model 
at applied loads of 24.1 and 29.4 kN. 

The axial strains for the two gages shown in Figure 
15 are shown in Figure 18. The strains and the failure 
load predicted by the model correlate well with the test 
results. However, the test results after the debond load 
are misleading because the load was immediately 
reversed after failure by the test operator, even if the 
specimen could carry additional loads. 
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Figure 18. Strains at two gage locations for test and 
models with and without ply damage 

The extensometer measurements for tests C4, C6, 
C8 and CIO conducted by Cvitkovich 23 and the 
predicted values are shown in Figure 19. The initiation 
of delamination is marked by the sharp breaks in the 
extensometer readings. After failure, the extensometer 
measurements are unreliable, as can be observed by 
comparing the different tests. The predicted response 
correlates well with the experimental results. The 
predicted debond load is just outside the range for all 
four tests. 



Figure 19. Predicted and experimental extensometer 
measurements. 


Concluding Remarks 

A method for the simulation of progressive 
delamination based on decohesion elements was 
presented. Decohesion elements are placed between 
layers of solid elements and they open in response to the 
loading situation. The onset of damage and the growth 


of delamination can be simulated without previous 
knowledge about the location, the size, and the direction 
of propagation of the delaminations. A strain softening 
law for mixed-mode delamination was proposed. The 
criterion uses a single state variable, the maximum 
relative displacement 6 ” ll “ , to track the damage at the 
interface under general loading conditions. The material 
properties required to define the element constitutive 
equation are the interlaminar fracture toughnesses and 
the corresponding strengths. 

Four examples were presented that test the 
accuracy of the method. Simulations of the DCB and 
ENF test represent cases of single-mode delamination. 
The MMB test that was simulated has equal Mode I and 
Mode II. The skin/stiffener debond problem was used 
as a test of the capabilities of the method in an 
important structural problem that exhibits non self- 
similar delamination growth. The examples analyzed 
indicate that the method of decohesion elements can 
predict the strength of composite structures that exhibit 
progressive delamination. 
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